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The a(4) Scheme — A High Order Neutrally Stable CESE Solver 


Sin-Chung Chang 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 

Abstract 

The CESE development is driven by a belief that a solver should (i) enforce conservation 
laws in both space and time, and (ii) be built from a non-dissipative (i.e., neutrally stable) 
core scheme so that the numerical dissipation can be controlled effectively. To provide a solid 
foundation for a systematic CESE development of high order schemes, in this paper we describe 
a new high order (4- 5th order) and neutrally stable CESE solver of the advection equation 
du/dt + adu/dx = 0. The space-time stencil of this two-level explicit scheme is formed by one 
point at the upper time level and two points at the lower time level. Because it is associated 
with four independent mesh variables u (w x )™, (u xx )^ and (w X xx)j (the numerical analogues 
of u, du/dx , d 2 u/dx 2 , and d 3 u/dx 3 , respectively) and four equations per mesh point, the new 
scheme is referred to as the a( 4) scheme. As in the case of other similar CESE neutrally stable 
solvers, the a (4) scheme enforces conservation laws in space-time locally and globally, and it has 

the basic, forward marching, and backward marching forms. Assuming \v\ ^ 3 (y = f aAt/Ax ), 
these forms are equivalent and satisfy a space-time inversion (STI) invariant property which 
is shared by the advection equation. Based on the concept of STI invariance, a set of algebraic 
relations is developed and used to prove that the a (4) scheme must be neutrally stable when 
it is stable. Numerically, it has been established that the scheme is stable if \v\ < 1/3. 

1. Introduction 

The space-time conservation element and solution element (CESE) method is a high- resolution and 
genuinely multidimensional method for solving conservation laws [1-73]. Its nontraditional features include: 
(i) a unified treatment of space and time; (ii) the introduction of conservation elements (CEs) and solution 
elements (SEs) as the vehicles for enforcing space-time flux conservation; (iii) a novel time marching strategy 
that has a space-time staggered stencil at its core and, as such, fluxes at an interface can be evaluated 
without using any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability 
to capture shocks without using Riemann solvers); (iv) the requirement that each scheme be built from a 
non-dissipative core scheme and, as a result, the numerical dissipation can be controlled effectively; and (v) 
the fact that mesh values of the physical dependent variables and their spatial derivatives are considered as 
independent marching variables to be solve for simultaneously. 

Without using flux- splitting or other special techniques, since its inception in 1991 [1] the unstructured- 
mesh compatible CESE method has been used to obtain numerous accurate ID, 2D and 3D steady and 
unsteady flow solutions with Mach numbers ranging from 0.0028 to 10 [51]. The physical phenomena 
modeled include traveling and interacting shocks, acoustic waves, vortex shedding, viscous flows, detonation 
waves, cavitation, flows in fluid film bearings, heat conduction with melting and/or freezing, electrodynamics, 
MHD vortex, hydraulic jump, crystal growth, and chromatographic problems [3-73]. In particular, the rather 
unique capability of the CESE method to resolve both strong shocks and small disturbances (e.g., acoustic 
waves) simultaneously [13,15,16] makes it an effective tool for attacking computational aeroacoustics (CAA) 
problems. Note that the fact that second-order CESE schemes can solve CAA problems accurately is an 
exception to the commonly- held belief that a second-order scheme is not adequate for solving CAA problems. 
Also note that, while numerical dissipation is needed for shock capturing, it may also result in annihilation of 
small disturbances. Thus a solver that can handle both strong shocks and small disturbances simultaneously 
must be able to overcome this difficulty. 
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In spite of its nontraditional features and potent capabilities, the core ideas of the CESE method are 
simple. In fact, all of its key features are the inescapable results of an honest pursuit driven by these 
simple ideas. The first and foremost is the belief that the method must be solid in physics. As such, in 
the CESE development, conservation laws are enforced locally and globally in their natural space-time unity 
forms for ID, 2D and 3D cases. Moreover, because direct physical interaction generally occurs only among 
the immediate neighbors, use of the simplest stencil also becomes a CESE requirement. Obviously, this 
requirement is also very helpful in simplifying boundary-condition implementation. 

The second idea is derived from the realization that stability and accuracy are two competing issues 
in time- accurate computations, i.e., too much numerical dissipation would degrade accuracy while too little 
of it will cause instability. In other words, to meet both accuracy and stability requirements, computation 
must be performed away from the edge (“cliff”) of instability but not too far from it. This represents a real 
dilemma in numerical method development. As an example, schemes with high-order accuracy generally has 
high accuracy and low numerical dissipation. However, it is susceptible to instability. In fact, in dealing with 
complicated real-world problems, stability of these schemes often is difficult to maintain without resorting 
to ad hoc treatments. To confront this issue head-on, in CESE development, it is required that a solver 
be built from a non- dissipative (i.e., neutrally stable) core scheme. By definition, computations involving 
a neutrally stable scheme are performed right on the edge of instability and therefore the numerical results 
generated are non-dissipative. As such numerical dissipation can be controlled effectively if the deviation of 
a solver from its non-dissipative core scheme can be adjusted using some built-in parameters. Note that the 
above idea also plays an essential role in the recent successful development of a family of Courant number 
insensitive schemes [59,61,64,65,67]. 

Other CESE ideas are: (i) the flux at an interface be evaluated in a simple and consistent manner; 
(ii) genuinely multidimensional schemes be built as simple, consistent and straightforward extensions of 
ID schemes; (iii) triangular and tetrahedral meshes be used in 2D and 3D cases, respectively, so that the 
method is compatible to the simplest unstructured meshes and thus can be used to solve problems with 
complex geometries; and (iv) logical structures and approximation techniques used be as simple as possible, 
and special techniques that has only limited applicability and may cause undesirable side effects be avoided. 
Fortunately for the CESE development, as it turns out, the realization of the above lesser ideas (i)-(iv) 
follows effortlessly from that of the first two core ideas. 

The first model equation considered in the CESE development is the simple convection equation 


du du 


(i.i) 


where the advection speed a 7 ^ 0 is a constant. Let x\ = x, and 
of a two-dimensional Euclidean space E<i. Then, because Eq. (1.1) 

h d = (au,u), Gauss’ divergence theorem in the space-time E<i implies 
of the integral conservation law 


/ 

JS(V ) 


h • ds = 0 


= t be considered as the coordinates 
can be expressed as V • h = 0 with 

that Eq. (1.1) is the differential form 

( 1 . 2 ) 


As depicted in Fig. 1 , here (i) S(y) is the boundary of an arbitrary space-time region V in ^ 2 , and (ii) 
ds = da n with da and ft, respectively, being the area and the unit outward normal vector of a surface element 
on S(V). Note that: (i) because h ■ ds is the space-time flux of h leaving the region V through the surface 
element ds, Eq. (1.2) simply states that the total space-time flux of h leaving V through *S(I/) vanishes; 
(ii) in E 2 , da is the length of a line segment on the simple closed curve S'(V^); and (iii) all mathematical 
operations can be carried out as though E<i were an ordinary two-dimensional Euclidean space. 

It is well known that a solution to Eq. ( 1 . 1 ) represents non-dissipative data propagation along its 
characteristic lines defined by dx/dt = a. Moreover, Eq. ( 1 . 1 ) is invariant under space-time inversion (STI), 
i.e., it transforms back to itself if x and t are replaced by —x and —t, respectively. (In physics, STI invariance 
generally is referred to as PT invariance where P denotes a mirror- image or spatial- reflect ion operation while 
T denotes a time- reversal operation). Thus a solution to Eq. ( 1 . 1 ) possesses the following properties: (i) it 
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is completely determined by the data specified at an initial time level; (ii) its value at a space-time point 
has a finite domain of dependence (a point) at the initial time level; and (iii) the space-time inversion image 
of a solution to Eq. (1.1) is also a solution and vice versa. As such, in the initial CESE development, the 
focus is on the construction of an ideal core solver of Eq. (1.1) that enforces the conservation law Eq. (1.2) 
and also possesses all other properties of Eq. (1.1), i.e., it is a two-level, explicit, non-dissipative, and STI 
invariant solver. An in-depth account of this development and the resulting “a” scheme is given in [71]. As 
it turns out, the 2nd-order accurate a scheme (i) has a space-time stencil formed by one mesh point at the 
upper time level and two mesh points at the lower time level; and (ii) it is neutrally stable if v 2 < 1 where 
v = aAt/Ax. Also, at each space-time mesh point (j, n) , the a scheme is associated with two independent 
mesh variables u™ and ( u x (the numerical analogues of u and dujdx , respectively) and two equations. 

Until recently, with one exception (a three-level and 3rd-order accurate scheme reported on p. 80 of 
[1]), all CESE solvers of Eq. (1.1) are two- level and 2nd-order accurate extensions of the a scheme. To 
initiate a systematic development of CESE schemes with high order accuracy, two new high order accurate, 
conservation-law enforcing, and neutrally stable CESE solvers of Eq. (1.1) has been developed recently. 
Both solvers are explicit and involving two time levels. The space-time stencil of one of them is formed by 
one point at the upper time level and three points at the lower time level. Because it is associated with 
three independent mesh variables u™, (u x )™ and (u xx )™ (the numerical analogues of u, dujdx , and d 2 u/dx 2 , 
respectively) and three equations per mesh point, the scheme is referred to as the a(3) scheme in [72]. On the 
other hand, the space-time stencil of the second scheme to be described here is formed by one point at the 
upper time level and only two points at the lower time level. Because it is associated with four independent 
mesh variables u™, ( u x )™, (u xx )‘j , and (w xxx )^ (the numerical analogues of u, dujdx, d 2 u/dx 2 , and d 3 u/dx 3 , 
respectively) and four equations per mesh point, hereafter the new scheme is referred to as the a (4) scheme. 

2. The a( 4) scheme 

To proceed, consider the set of space-time staggered mesh points (j,n) (dots in Fig. 2(a)), where 

Qi = f {(j, n)\j, n = 0, ±1, ±2, d=3, . . . , and (j + n) is an odd integer} (2.1) 

Each (j, n) £ fti is associated with a solution element, i.e., SE(j, n) (see Fig. 2(b)). Let points 77, G, L, 
and M (marked by small open circles in Figs. 2(b)-2(f)) be the midpoints of the line segments AF, AG, 
A/, and ATT, respectively. Then, by definition, SE(j, n) is the interior of the space-time region bounded by 
a dashed curve depicted in Fig. 2(b). It includes a horizontal line segment ED, a vertical line segment BJ , 
two inclining line segments HL and GM, and their immediate neighborhood. 

At this juncture, the reader is warned that the notation used here may differ from that used in previous 
CESE papers [1-70]. In particular, (i) the mesh indices j and n are only allowed to be whole integers here; 
and (ii) the spatial and temporal intervals that are denoted by Axj 2 and At/2, respectively, in [1-70] are 
denoted by ax and At, respectively, here. These changes are introduced so that the a (4) scheme can be 
compared with the a (3) scheme on the same footing. Note that the a (3) scheme, like most established 
schemes, is constructed over a set of mesh points which are not staggered in space-time. 

Let ( x,t ) £ SE(j, n). Then Eqs. (1.1) and (1.2) will be simulated numerically assuming that u(x,t) and 
h(x,t), respectively, are approximated by 

u*(x,t;j,n ) = f u"+ (u x )y(x-Xj) + (u t )"(t - t") + ^{u xx )](x - xj ) 2 + (u xt )j(x - Xj)(t - t n ) 

+ - t n f + !(«***)?(* - Xj ) 3 + \{u xxt y]{x - Xj) 2 (t - t n ) (2.2) 

+ ^(Uxtt)j(x - Xj)(t ~ t n ) 2 + ^(u ttt )](t - t n ) 3 

and 

h*(x,t;j,n ) = f ( au*(x,t;j,n ), u*(x,t;j,n)) (2.3) 
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Note that: (i) u £, ( u x )?, (u t )?, (u M )J, (w x t)j , (u«)” (u IM )J, (w xx t)j » and (wttt)j are constants 

in SE(j, n), and the numerical analogues of the values of w, du/dx , du/dt , d 2 u/dx 2 , d 2 u/dxdt , d 2 u/dt 2 , 
d 3 w/<9z 3 , d 3 u/dx 2 dt , d 3 u/dxdt 2 , and d 3 u/dt 3 at the mesh point (j, n), respectively; (ii) (a^i™) are the 
coordinates of the mesh point (j, n) where 2Cj = jA# and i n = nAt ; (iii) u*(x,t;j,n) represents a 3rd-order 
Taylor’s approximation of u; and (iv) Eq. (2.3) is the numerical analogy of the definition h = (au,u). 

For any ( j, n) £ ffy, let w = w*(:r,t ; j, n) satisfy Eq. (1.1) for all (x,t) £ SE(j, n). Then one has 


= -a(u x )J, (u xt )j = -a(u xx )j, (utt)j = a 2 (u xx )", 

xxt)j = &(^XXx)j ) ( 'U'xttjj = ® (^xxx)j ) = ® (^xxx)j 


(2.4) 


Substituting Eq. (2.4) into Eq. (2.2), one has 


u*(x,t-,j,n) = uj+ {u x yy [(x-Xj) - a(t-t n )] + huxx)" [(z - Zj) - a(t-t ")] 2 


+ h«xxx)" [(a: ~ £j) -a(t- t n )] 3 


(j,n)G^i (2.5) 


i.e., Uj , (w x )^, (w X x)j 5 and (w xxx )j are the only independent mesh variables associated with (j, n). 

With the above preliminaries, next we describe the basic form of the a (4) scheme. 

2.1. The basic form of the a (4) scheme 

Let Ei be divided into non- overlapping space-time triangular regions (see Fig. 2(a)) referred to as 
conservation elements (CEs). As depicted in Figs. 2(c)-2(f), (i) each ( j, n) £ f2i is assigned with four CEs, 
i.e., CE(j, n;^), £ = 1,2, 3, 4; (ii) each CE represents a right triangle with the end points of its hypotenuse 
£ ffy but not the third vertex; and (iii) the space-time E 2 can be filled by CE (j,n;£), £ = 1,2, 3, 4; and 
(j, n) £ The a( 4) scheme will be constructed by assuming that the flux of h* conserves over CEs, i.e., 


/ 

JS(CE(j,n ,£ )) 


h* ■ ds = 0, 


^=1,2, 3, 4; (j,n) £fi 1 


(2.6) 


Using the special case £ = 1 as an example, how Eq. (2.6) can be turned into a set of relations linking 
the mesh variables at two diagonally opposite neighboring mesh points will be explained step-by-step in the 
following remarks: 

(a) By definition, on S(CE(j , n; 1)) (i.e., the boundary of CE (j, n; 1)), the line segments AD and AG belong 
to SE(j, n) while CD and CG belong to SE(j + 1, n — 1). Note that, strictly speaking, points D and G 
do not belong to either SE(j, n) or SE (j + l,n — 1). This fact, however, does not pose a problem for 
flux evaluation over S(CE(j , n; 1)) because the values of h* at isolated points do not contribute to the 
flux of h* over a finite line segment. 

(b) The straight line passing through points A and C can be defined by 

t = t n — — (x — Xj) or t = t n ~ 1 — —(x — Xj+i) (2.7) 

AX AX 


(c) For CE(j, n; 1), the outward unit normal vectors ft on AD, CD, and AC are 


ni = f (0, 1), n 2 d =(l,0), 


and 


def 

n 3 = - 


(a t, ax) 


V(a*) 2 + (ax) 2 


( 2 . 8 ) 


respectively. 

(d) Obviously, (i) the length of the line segment joining any two neighboring points (x,t n ) and (x + dx , t n ) 
on AD is dcr = \dx\\ and (ii) the length of the line segment joining any two neighboring points (xj+i,£) 
and (fljj*_|_i,t + dt) on CD is da = \dt\. 
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On the other hand, because \dt/dx\ = At/ ax for any two neighboring points (x,t) and (x A dx,t + dt) 
on AC, the length of the line segment joining these two points is 

da = vW + W = V(At) 2 + (az) 2 ■ (\dx\/Ax) (2.9) 

(e) With the aid of (i) Eq. (2.3), (ii) the preliminaries given in the above items (c) and (d), and (iii) the 
relation ds = daft, one concludes that, for CE (j, n; 1), 

{ u*\dx\ on AD 

au*\dt\ on CD (2.10) 

— (1 + v)u*\dx\ on AC 


Here v = f aAt/ ax is the Courant number. 

(f) With the aid of Eq. (2.10) and the comments made in the above items (a) and (b), one can see that, 
for the case t = 1, Eq. (2.6) <=$ 


pXj+i rt 

/ u*(x,t n ; j, n)dx A / au*(xj+i,t; j + 1, n — l)dt 
J Xj 

r x J+ V 2 /a t \ 

— (1 + v) / u*[x,t n (x — Xj); j, n)dx (j,n)eQi (2.11) 

J Xj \ ax J 

f Xj+1 / At \ 

— (1 + v) / (x — Xj+i);j A l,n — 1 )dx = 0 

1/. V A:r J 


9 x j+l(2 

where (i) the symbol “^=>” is used as a shorthand for the statement “if and only if”, and (ii) £j+i /2 = f 
Xj A (ax/2). 

(g) Let 

(«*)? = f ^ (uxx)] = and (u^)J = ^/~(u xxx )?, y,n)€fti (2.12) 

Then, with the aid of Eq. (2.12), the expression obtained by substituting Eq. (2.5) into (2.11) can be 
cast into the following form: 


,, , f ^ + 3 v 2 + 4v + 7 , (y + Z)(v 1 + 2v + 5) i" 

(1 Vj \ U A 2 'U'x A g U xx A n A 'Llxxx 


24 


J 3 


{l-v) 


3v A 1 7z/ 2 + 4zy + l (3v A l)(5i^ 2 + 2v + 1) i n ~ 1 

W ~ Ux H ^ %x %xx 


(j,n) Gfii (2.13) 


6 


24 


3 + 1 


To simplify notation, in the above and hereafter we adopt a convention that can be explained using an 
expression on the left side of Eq. (2.13) as an example, i.e., 


z^ + 3 i/ 2 + 4i' + 7 (v A 3)fy 2 + 2v + 5) ~\ n 

U H n U X H ^ U XX H fTi u xxx 

2 o 24 J j 

= «J + + l ' 2+ ^ + 7 (. M )? + (- + 3)(^+2- + 5) (um8 )„ 


24 
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Similarly, for the cases l = 2,3,4, Eq. (2.6) implies that 


= (!-*') 


3^+1 7i' 2 +4z; + l (3u+ 1)(5i ^ 2 + 2i/ + 1) i» 

W “I" ~ “I" 7 lAxxx 

Z O 24 J j 


v + 3 v 2 + 4u + 7 (v + 3)(v 2 +2i/ + 5) "| 71-1 

U — U x ”1" 7 ^xx J77 Uxxx 


(j, n) G Hi 


(2.14) 


24 


i+i 


. r i/ — 3 z/ 2 -4z/+7 (i/ - 3)(i^ 2 - 2v + 5) i n 

(1 “t“ ^) \^3> “t“ ^ H g “I" ^ %xx 

. r 3^ — 1 7 i/ 2 -4z/+ 1 (3^ — l)(5i^ 2 — 2i/ + 1) i"- 1 

(1 “I - 2 4 g %x ^ 'U'xxx _ ^ 


(j»€ni (2.15) 


and 


(! + ")[ 


3i/-l 7^ 2 — 4i/ + 1 (3i^ — l)(5i^ 2 — 2z/ + 1) 

W 4 » 4 7 4 pTx 'Uxxx 


(l+*0 


1^-3 i/ 2 - 4i/ 4- 7 

^ Pi ^X 4 7 ^XI 


24 


(^ — 3) (i^ 2 — 2i^ 4- 5) -in- 1 

^xxx 


(j, n) G Hi 


(2.16) 


24 


j~ 1 


respectively. 

At this juncture, note that: 

(a) Because 

<9^ ax <9-u <9 2 w (ax ) 2 d 2 u d 3 u (ax)3 d 3 u _ def x 

dx 2 <9x 5 dx 2 4 dx 2 5 an dx 3 8 <9x 3 1 X ax/2 

the normalized parameters (tlx)™, (^xx)jN and (w^xx)™ can be interpreted as the numerical analogues 
of the values at (j, n) of the first, second, and third derivatives of u with respect to the normalized 
coordinate x. 

(b) Note that: (i) the vector h* at any horizontal or vertical interface separating two neighboring CEs is 
evaluated using the information from a single SE; (ii) the vector h* at one half of any inclining interface 
separating two neighboring CEs is evaluated using the information from a single SE while that at another 
half is evaluated using the information from another SE; and (iii) the unit outward normal vector on the 
surface element pointing outward from one of two neighboring CEs sharing the same element is exactly 
the negative of that pointing outward from another CE. Thus one concludes that the flux leaving one 
of the two neighboring CEs through the interface they share is the negative of that leaving another 
CE through the same interface. Due to this interface flux cancelation and the fact that the CEs are 
nonoverlapping and can fill the space-time E 2 , the local conservation relations Eq. (2.6) lead to a global 
conservation relation, i.e., the total flux of h* leaving the boundary of any space-time region that is the 
union of any combination of CEs vanishes. 


Let 1 — v 7^ 0 and 1 4- v ^ 0, i.e., 

i/ 2 ^ 1 (2.17) 

Then Eqs. (2.13)-(2.16) can be simplified by eliminating (i) the common factors (1 — v) on the both sides of 
each of Eqs. (2.13) and (2.14); and (ii) the common factors (1 4- z^) on the both sides of each of Eqs. (2.15) 
and (2.16). By adding the simplified forms of Eqs. (2.13) and (2.14) together and then subtracting one of 
them from another, one has 


[u + (1 + v) 


, 2(l + zy + z/ 2 )^ ( (l + i/)(l + i/ 2 ) 1" 

u x 4“ u xx 4" u xxx 


M , .A.. I 2(l+I/ + J/ 2 )_ ( 1 + I ')( 1 +*' 2 )„. I”" 1 

U (1 + V)U X + U xx V'xxx 


(j»6fli (2.18) 


3 + 1 
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and 


„ . 7z/ 2 + 10z/ + 7 1" 

U X + (1 + l^j'U'XX "t" -y n 'U'XXX | 


12 




( v 7z/ 2 + 101/ + 7 

(1 H - )'LL XX H ^2 'LLxxx 


n— 1 

j 

J+l 


(j,n)Gfii (2.19) 


Note that Eq. (2.19) has been further simplified by eliminating the common factors (1 — v) which appear 
on the both sides of the original subtraction result. Similarly, the simplified forms of Eqs. (2.15) and (2.16) 
imply that 


[“-(!- v ) 


2(l-i/ + z/ 2 ) (1 - i/)(l + z/ 2 ) -|« 

'LL X "I - n 'LL xx I 'LL XXX 

o O -I j 


,, . 2(1-1/ + 1/ 2 ) ( 1 -z/)(l + z/ 2 ) 1"-1 

'LL + (1 L/jU x “h _ Wjjx + „ | 


(j, njefil 


( 2 . 20 ) 


J J-l 


and 


(l ^)^xx “h 


7z/ 2 - 10i/ + 7 
12 1 


“I" (l lp^)u xx ~\~ 


7v 2 - 10i/ + 7 
12 ~~ 


^xxx 


n— 1 

j 

i-i 


(j,n)€fl i (2.21) 


Note that Eqs. (2.18) and (2.19) are stricter than Eqs. (2.13) and (2.14) in the sense that the former imply 
the latter for any v but the latter imply the former only if an extra condition (i.e., 1 — v ^ 0 for this case) is 
imposed. Similarly, Eqs. (2.20) and (2.21) are stricter than Eqs. (2.15) and (2.16). Hereafter, by definition, 
the system of equations formed by Eqs. (2. 18) -(2. 21) is referred to as the basic form of the a (4) scheme. In 
the following, we derive the forward marching form. 

2.2. The forward marching form of the a (4) scheme 

To simply notation, temporally the expressions on the right sides of Eqs. (2. 18) -(2. 21) will be denoted 
by s i, S2, S3, and S4, respectively. Then, by adding Eqs. (2.18) and (2.20) together and then subtracting 
one of them from another, one has 


2(1 + z/ 2 ) i/(l + z/ 2 ) 

'LL ”|” V'LL X "T ” %x n *LLxxx 


= 2( s i + s 3 )> ( j , n)en 1 


and 


2v 1 + v 2 

'L L x + ~~q~'LL xx “h ~ U xxx 


= 2 (^1 33), (j,n)€Sl 1 


In turn, by subtracting Eq. (2.23) from both Eqs. (2.19) and (2.21), one has 


v -|- 3 


3z/ + l 

^xx H“ A 'LL xxx 


Si- Sz), 1 


and 


1 / — 3 


'LLxx T 


3z/-l 


'LLxxx 


J 3 


S4-2( s i _s 3). 1 


Let v — 3^0 and v + 3 ^ 0, i.e., 
Then Eqs. (2.24) and (2.25) 


v 2 ^ 9 


~, n _ 6 s 2 6 s 4 1 18(si - S3) 

['LLxxx) j — | q o' 9 n ’ 

J V + 3 V — 3 V A — 9 


(j, n) G Oi 


( 2 . 22 ) 

(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 
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( \n _ 3(3i/-1)«2 3(3z/ + 1)s4 15z/(si - S3) 

{U ™ )j ~ 2(i/ + 3) + 2(i/ -3) i/ 2 -9 ’ 

Next, by substituting Eqs. (2.27) and (2.28) into Eq. (2.23), one has 

_ (z/ + 1)(z/-2)s 2 _ {v- 1)(i/ + 2)s 4 3(3i/ 2 - 7)(si - S3) 

v + 3 i/ - 3 + 2(i/ 2 - 9) 

Finally, by substituting Eqs. (2.27)-(2.29) into Eq. (2.22), one arrives at 


C?\ ^0 ^ ^1 


(j, n) e 


,n_( 3l/ - 1 ) s 2 (3i/+1)s 4 1 [ I/(l/2 — 29)1 1 r i/(i/2-29)1 _ ~ /o on', 

^““T+3 ^3~ + 2 1 ^9~ Sl + 2 1+ ~^9~ S3 ’ (j>n)eni (2 - 30) 


iv 2 — 9 


Hereafter , in Eqs. (2.27)-(2.30) and all other equations derived from them, Eq. (2.26) will be assumed 
implicitly. 

To proceed, let 


n -\ ( u *)j 

q{3 ' ] («xx)? 

V (Uxxx ) j 


(j, n) G ft 1 


d ^ f ( Sl ) and s_ d ^ f ( S3 


C+{v) 


def / 1 — (1 + Z/) (2/3) (1 + V + Z/2) —(1 + Z/)(l + Z/ 2 )/3 

“ lo 1 -(1 + z/) (7z/ 2 + 10z/ + 7)/12 


C_(z/) = 


def / 1 1 — z/ (2/3) (1 — z/ + z/ 2 ) (1 — z/)(l + z/ 2 )/3 


0 1 


(7z/ 2 - 10z/ + 7)/12 


D + (z/) S 


/{l - - 29)/(^ 2 - 9)]}/2 (3z/ - l)/(z/ + 3) 

def (3/2)(3z/ 2 - 7)/(z/ 2 — 9) (z/ + l)(z/ — 2)/(z/ + 3) 

— 15z//(z/ 2 — 9) — (3/2)(3z/— l)/(z/ + 3) 

V 18/ (z/ 2 — 9) 6/(z/ + 3) 


£>_(z/) = 


(1 + K^ 2 - 29)/ (z/ 2 - 9)]}/2 — (3z/ + l)/(z/ - 3) 

def — (3/2)(3z/ 2 — 7)/(z/ 2 — 9) — (z/ — l)(z/ + 2)/(z/ — 3) 

15z//(z/ 2 — 9) (3/2)(3z/+l)/(z/-3) 

V — 18/(z/ 2 - 9) — 6/ (z/ - 3) 


Then: (i) by the definition of si, s 2 , S3, and s 4 , one has 

s+ = C+ (z/) $0 + 1, n - 1) and s_ = C_ (z/) (ft; - 1, n - 1), ( j , n) G 


(2.36) 


and (ii) Eqs. (2.27)-(2.30) can be cast into the form 

q(j, n) = D + (z/) s+ + D_ ( 1 /) s_ , 


(j, n) G fli 


By substituting Eqs. (2.37) into (2.38), Eqs. (2.27)-(2.30) can be cast into the matrix form 


q{j, n) = Q + ( v) q(j + 1, n - 1) + Q-(v) q{j - 1, n - 1), 


0, n) G Cli 


Q+( v) = f D. |_(z/) C+{ v) and Q-( v) = f D_(z/) C_(z/) 
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Hereafter the system of equations formed by Eqs. (2.27)-(2.30) or its equivalent Eq. (2.39) is referred to as 
the forward marching form of the a( 4) scheme. According to the above derivation, the forward marching 
form is a result of the basic form and the assumption Eq. (2.26). Also note that here Eq. (2.40) is cast into 
a form different from the forward marching form of the a scheme (e.g., Eq. (3.48) in [71]), i.e., the symbols 
Q+(z') and Q-(v) used in the a scheme are replaced in the a( 4) scheme by Q- (v) and Q+( v), respectively. 

As a preliminary to a later development, next we will take a side tour and introduce the concept of 
invariance under space-time inversion. 

2.3. Invariance under space-time inversion 

Let u = u(x,t) be a solution to Eq. (1.1) in the domain — oo < x, t < +oo, i.e., 


Let 

and 

Then (i) Eq. (2.41) <^> 


du(x,t) du(x,t) _ n 
— dt~ +a — 9 ^=°’ 


— oo < x,t < +oo 


/ def 
X = 


—X 


and t' = -t 


u(x,t) = u(—x,—t) 


(2.41) 

(2.42) 

(2.43) 


du(x',t') duix' ,t') 

-M^ + a -8*- 30 ’ 


— oo < x f ,t' < -f oo 


and (ii) 

d___d_ _d_ _ _d_ 

dt' ~ 9t an dx' ~ dx 

Thus one concludes that Eq. (2.41) 


du(x.t) du(x,t) 

-g^ + a S^ s0 ' 


—oo < X, t < -t-oo 


(2.44) 


(2.45) 


(2.46) 


In other words, if u = u(x,t) is a solution to Eq. (1.1), so must be u = u(x,t) and vice versa. Because the 
one-to-one mapping 

(x, t) <-► (— x, — t), — oo < x,t < +oo (2.47) 

represents a space-time inversion (STI) operation, hereafter (i) a pair of functions such as u and u will be 
referred to as the STI images of each other; and (ii) a partial differential equation (PDE) such as Eq. (1.1) 
is said to be STI invariant if the STI image of a solution is also a solution and vice versa. 

Next let 


,(M) 


(x, t) 


def d k+£ u(x,t ) 


and u^ k,t ^{x,t) 


def d k+e u(x,t) 


dx k 0t* ~ 8x k dt i ’ 

Then, with the aid of the chain rule, Eqs. (2.42), (2.43), and (2.48) imply that 


oo < x,t < H-oo; k,£ = 0,1,2,... 

(2.48) 


M (x t) = g fc+< «(-a:, -t) = k+e dyMyX) 

K,) dx k dt e K ’ dx ,k dt ,e — oo < x,t < +oo; = 0, 1, 2, . . . (2.49) 

= (-1 ) k+e u ik ’Q(x',t') = (-1 ) k+e u( k ’ e \-x,-t) 


i.e., 


u( k ’ e \x,t) 


( u( k,e \—x, — t ) if ( k + i) is even 
\ —u^’^ (—x, — t ) if ( k + t) is odd 


(2.50) 
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According to Eq. (2.48), = u and = u. Thus Eq. (2.43) is a special case of Eq. (2.49) with 

k = £ = 0. 

In the following, the concept of STI invariance will be introduced for the a( 4) scheme. As a preliminary, 
note that: (i) (j, n) £ fti (— j, — n) £ fii; (ii) 

0\ n) <-► (-j,-n) (2.51) 

is the numerical analogue of the STI mapping Eq. (2.47); and (iii) , (w x )^ , (w*)™, (w xx )^, (uxt)j , (^tt)j 
(^xxx)j ? (^xxt)j ? (uxtt )™ , and (^m)j are the numerical analogues of the values of w, du/dx , du/dt , d 2 u/dx 2 , 
d 2 u/dxdt , d 2 u/dt 2 , d 3 u/dx 3 , d 3 u/dx 2 dt , d 3 u/dxdt 2 , and d 3 u/dt 3 at the mesh point (j, n), respectively. 
Thus, motivated by Eq. (2.50), the one-to-one mapping 

u™ <-+ ul"; (w*)j <-► -(^x)Ij ; (w*)j «-► -(«*)!" ; (w**)? «-► fa**)!” 

faxt)j faxi)_j j faii)j fatt)_j 5 fa***)j faxxx)_j (j? ^ ^1 (2.52) 

faxxi)j <"♦ -fa**t)lj 5 fa*tt)j «"► -fa*tt)l”; (wm)j «"► -fattt)l” 
is taken as the numerical analogue of the one-to-one mapping 

u( fc ^)(x,t) <-► fi( fc, ^(x, i), — oo < a:, t < +oo; fc,^ = 0, 1,2,3 (2.53) 

For the independent mesh variables, by using Eq. (2.12), it is seen that the mapping Eq. (2.52) reduces 


/ Uj \ ( U -1 

(«.)? 1 „ -(«*)=; 
(^xx)" , (Wxi)_j 

V(^xii)"/ \—(Uxxx) Ij/ 

With the aid of Eq. (2.31), Eq. (2.54) can be expressed as 

0, n) G Hi 

(2.54) 

<f0» «-*■ Uq(-j,-n), 

0, n) G Hi 

(2.55) 

where 

/I 0 0 0 \ 

jj def 0 — 10 0 

^ ” 0 0 1 0 

Vo 0 0 -1/ 

The matrix U is unitary. In fact it is a real matrix with 

(2.56) 

II 

Cj 

1 


(2.57) 


Hereafter (i) M -1 denotes the inverse of any nonsingular square matrix M; (ii) for each (j, n), Uq[—j, —n) is 
referred to as the STI image of q(j, n); and (iii) the set formed by Uq(—j, — n), (j, n) £ is also referred to 
as the image of the set formed by q(j , n), (j, n) £ Sli. According to Eq. (2.57), q(j , n) = UUq{— (— j), — (— n)). 
Thus q(j, n) is the STI image of Uq(—j,—n) as an individual (j, n) or as the set defined over fii. In the 
following, we will show that the system of equations defined by each of Eqs. (2.18)-(2.21) is STI invariant, 
i.e., the system maps onto an equivalent system under the mapping Eq. (2.54). 

As an example, consider the system defined by Eq. (2.18). Under the mapping Eq. (2.54), Eq. (2.18) 
maps onto 


M , ..i„. , 2(1 -\-v + v 2 ) __ (1 + IV)(1 + ry 2 ) ■ 

U (1 “h v)u x T „ Wju „ 'U'xxx 


~3 


.. , „ , , 2(1 -\-v-\- v 2 ) _ , (l + i/)(l + i/ 2 )_ 

u “h (1 v)u x T „ T 'U'xxx 


-(n-l) 

-0+1) 


O’) n) G fil 


(2.58) 
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At this juncture, note that, in addition to changing the sign of each Ux and mapping Eq. (2.54) 

requires that the upper and lower indices j, n, j + 1, and n — 1 in Eq. (2.18) be replaced by their negatives, 
respectively. This is different from simply replacing the symbols j and n everywhere with —j and —n, 
respectively. Moreover, to simplify argument, hereafter system B is referred to as the STI image of system 
A if A maps onto B under the mapping Eq. (2.54), e.g., the system Eq. (2.58) is the STI image of Eq. (2.18). 
Let 

j* = f — j — 1 and n* d = —n + 1, (j, n) £ ffy (2.59) 

Then, by using the fact that (j* + n*) + (j + n) = 0 and therefore fy*,n*) € ffy (j, n) £ ffy, Eq. (2.58) 

can be cast into the form 


\ 2(1 + i/ + v 2 ) (1 + v){l + v 2 ) 

W + (1 + H - « 'LL xx "h q 'Llxxx 


„ , , 2(1 + + 1 ^ 2 ) {l + v){l + v 2 )~ V»*-i 

W (1 + L/JU X “h ^ 'LLxx q LLxxx 

O O 


0f\O Gfii (2.60) 


By comparing Eqs. (2.18) and (2.60), one can see that the system Eq. (2.18) is identical to its STI image 
Eq. (2.58) (which is identical to Eq. (2.60)). Thus, under the mapping Eq. (2.54), Eq. (2.18) maps onto 
itself, i.e., the system Eq. (2.18) is STI invariant. QED. 

The STI invariance of each of Eqs. (2.19)-(2.21) can be established in a similar manner. As such, the 
basic form of the a (4) scheme (which is formed by all component equations in Eqs. (2. 18) -(2.21)) is STI 
invariant. 

Let <f(j, n) = q Q (j , n), (j, n) £ Hi, be a solution to the basic form. Then, by substituting q(j , n) = q Q (j , n) 
into the basic form, one obtains a system of identities involving q Q (j, n), ( j, n) £ Due to the STI 
invariance of the basic form, the above system of identities is equivalent to that obtained by substituting 
= Uq 0 {— j,— n) into the basic form. As such q(j, n) = q Q (j, n), ( j, n) £ Sfy, represent a solution 
to the basic form q( j, n) = Uq 0 (— j, — n), ( j, n) £ ffy, represent another solution to the basic form. In 
other words, the STI image of a solution to the basic form is also a solution and vice versa. Obviously this 
conclusion is also valid for other STI invariant forms of the a (4) scheme. 

Next, we will establish the STI invariance of the forward marching form. As a preliminary, discussion 
of some basic concepts is in order. Note that for any set of variables X£, yi, £ = 1, 2, the conditions 


Xi + 2/1 = X 2 - 2/2 and x\ - yi = xi + yi 

(2.61) 

xi = xi and yi = -y 2 

Thus, the image of Eq. (2.61) under any one-to-one mapping 

(2.62) 

(xe,ye) <-*• i = 1,2 

i- e -, 

(2.63) 

+ 2/1 =x' 2 -y' 2 and x[ -y[=x 2 + y 2 
the image of Eq. (2.62) under the same mapping, i.e., 

(2.64) 

x[ = x 2 and y[ = -y 2 

(2.65) 


where the variables x f £ and y £, £ = 1,2, may or may not be related to l = 1,2. Moreover, in case 

that these two sets of variables are related, the condition Eq. (2.61) (or its equivalent Eq. (2.62)) may or 
may not be equivalent to the condition Eq. (2.64) (or its equivalent Eq. (2.65)). If the mapping Eq. (2.63) 
is such that Eq. (2.61) the image under this mapping (i.e., Eq. (2.64)), then Eq. (2.62) (the equivalent of 
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Eq. (2.61)) Eq. (2.65) (the equivalent of Eq. (2.64)). Eq. (2.63) with x\ = X£ and y £ = yt, £ = 1,2, is an 
example of such mapping while Eq. (2.63) with x\ = yt and y £ = X£, t = 1, 2, is not. 

The STI invariance of the forward marching form will be proved assuming Eq. (2.26) (the form is 
undefined if Eq. (2.26) is not valid). Note that: (i) the basic form of the a( 4) scheme <=> its forward marching 
form for any choice of q(j , n), (j, n) £ fii; and (ii) the STI images of the basic and forward marching forms, 
respectively, are obtained from the basic and forward marching forms through the mapping Eq. (2.55), i.e., 
through replacing q(j , n) in the basic form and the forward marching form with Uq(— j, — n), (j, n) £ Q. From 
the above observations and the illustration given in the last paragraph, one concludes that the STI image of 
the basic form that of the forward marching form. Because the basic form is STI invariant, i.e., the STI 
image of the basic form <=$ the basic form itself, Now we arrive at the conclusion that the forward marching 
form the basic form the STI image of the basic form <=$ the STI image of the forward marching form. 
Thus the forward marching form <=$ its STI image, i.e., the forward marching form is STI invariant. QED. 

2.4. The backward marching form of the a (4) scheme 

According to Eq. (2.55), the STI invariance of the forward marching form implies that Eq. (2.39) its 
STI image, i.e., 


Uq(-j,-n ) = Q+(v)Uq(-j - l,-n+ 1) + Q-{v)Uq{-j + 1 , -n + 1), (j,n) £ tti (2.66) 


By multiplying Eq. (2.66) from left using the matrix U and using Eq. (2.57), one concludes that Eq. (2.66) 
= Q-{v)q^-j- l,-n+l) + Q+{.v)^-j + l,-n+ 1), (2.67) 

where 

Q-{u) d = UQ+(u)U and Q+(v) = UQ-{v)U (2.68) 

By using Eqs. (2.33)-(2.36), (2.40), and (2.56), Eq. (2.68) implies that 

Q-(z') = D+(i/)C+(v) and Q+(^) = D~(v)C-(v) (2.69) 


where 


and 


C+(v) = C+(v)U 


(1 1 + u (2/3)(l + 1/ + v 2 ) (1 + u)(l + i/ 2 )/3 \ 

^0 -1 -(l + i/) -(7z/ 2 + 10u + 7)/12j 


C_( v) A =C-{v)U = 


b + {v)^UD + {v) = 


(1 — (1 — v) (2/3) (1 — v + v 2 ) -(l-^l + i/ 2 )^ \ 

VO -1 1-v ~{7u 2 - lQu+7)/12J 

/{I - ~ 29)/(^ 2 - 9)]}/2 (31/ - \)/{v + 3) X 

— (3/2) (3i^ 2 - 7)/ {y 2 - 9) -{v + l)(i/ - 2)/{y + 3) 

-15iV(i/ 2 - 9) — (3/2)(3i/ — l)/{v + 3) 

V — 18/(i/ 2 — 9) -6/(z/ + 3) / 


D-(v) = f UD_(v) 


(1 + Hu 2 - 29)/(i/ 2 - 9)]}/2 -(31/ + \)/{v - 3) \ 

(3/2) (3z/ 2 - 7)/{v 2 - 9) {v - l)(i/ + 2)/{v - 3) 
Yav/{y 2 — 9) (3/2)(3i/+1)/(i/-3) 

V 18/ (u 2 — 9) 6/(1 / -3) / 


(2.70) 

(2.71) 


(2.72) 


(2.73) 


Next, by replacing the “dummy” indices —j and — n everywhere in Eq. (2.67) with j and n, respectively, 
and using the fact that (— j,—n ) € fli & (j, n) e fii, one can see that the system Eq. (2.67) is identical to 
the system 


q{j,n) = Q+{v)q(j + l,n + l) + Q-(y)q(j- l,n + l), (j,n)eSl i (2.74) 


Because the mesh variables at (j, n) can be determined in terms of those at (j — 1, n + 1) and (j + 1, n + 1) 
using Eq. (2.74), hereafter Eq. (2.74) (which is equivalent to the forward marching forms of the a(4) scheme) 
will be referred to as the backward marching form of the a (4) scheme. 
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Eq. (2.74) was derived using the STI invariance of the forward marching form of the a( 4) scheme. 
Alternatively, it can also be derived from the basic form. Note that: (i) by replacing the “dummy” indices 
j and n everywhere in Eq. (2.18) with j — 1 and n + 1 and using the fact that (j, n) G fli (j + 1, n — 1) G 
fii (j — 1, n + 1) G fii, one can see that the system Eq. (2.18) is identical to the system 


u - (1 + v) 


Ux + 


Wtti 


U + (1 + v)Ux + 


2(l + i/ + i/ 2 ) (l + i/)(l + i/ 2 ) 

3 3 ’ 

2(l + i/+ i/ 2 ) (l + i/)(l + i/ 2 ) i«+ 1 


(j,n)Gfli (2.75) 


'U'XX 4~ 


^ULrn t t 


J J-l 


(ii) by replacing the indices j and n everywhere in Eq. (2.19) with j — 1 and n + 1 and using the fact that 
(j, n) GOi (j + 1, n — 1) (j - l,n+l) G 9 i, one can see that the system Eq. (2.19) is identical 

to the system 


'U'x ( 1 4“ ^) 'U'xx 4- 


7^ 2 + 10i^ + 7 I™ 


12 


WtT7 


4" (1 + v)llxx 4" 


7i^ 2 + 10v + 7 
12 ~ 


'U'XXX 


- 71+1 

5 

-I j-1 


0*,n)Gfii (2.76) 


(iii) by replacing the indices j and n everywhere in Eq. (2.20) with j + 1 and n + 1 and using the fact that 
(j, n) G 9i ^ (j — l,7i — 1) G fii -O- (j + 1, n + 1) G fii, one can see that the system Eq. (2.20) is identical 
to the system 


. 2(1 -v + v 2 ) (1 - iA(l + i/ 2 ) 

W 4” (1 l^j'U'x 4” q 'U'xx 4” 7) 'U'xxx 

o O J j 

2(1-1/ + 1/ 2 ) (l-i/)(l + i/ 2 ) 1»»+1 

U (1 V)U X + Wjj Uxxx 

3 3 Jj+i 


(j,n) G 


(2.77) 


and (iv) by replacing the indices j and n everywhere in Eq. (2.21) with j + 1 and n + 1 and using the fact 
that (j, n) G 9i ^ (j - l,n — 1) G (j + l,n 4- 1) £ ^i, one can see that the system Eq. (2.21) is 

identical to the system 


'tlx 4- (1 v^iixx 4“ 


7v 2 - 10i^ + 7 
12~” 


-'U'xxx 


J 3 


U, x (1 v'j'llxx 4” 


7z/ 2 - 10z/ + 7 
12 ~ 


Uxxx 


71+1 

-I j + 1 ’ 


(j,n)Gfi 1 (2.78) 


As such Eqs. (2.75)-(2.78) are equivalent to Eqs. (2.18)-(2.21), respectively. 

For each (j, n) Gfti, Eqs. (2.18)-(2.21) form a linear system of four equations for the four mesh variables 
u'j, (^x)j 5 (^xx)j ? and (^xxx)j ■ Eqs. (2. 75) -(2. 78) form another system. Moreover, one can see that, under 
the mesh variable mapping 


q(j , n) Uq(j, n), <f(j + 1, n - 1) Uq(j - 1, n + 1) and <f(j - l,n - 1) <-► g(j + 1, n + 1) (2.79) 

Eqs. (2.75)-(2.78), respectively, are the images of Eqs. (2.18)-(2.21) and vice versa. By using the concept 
introduced earlier in a discussion involving Eqs. (2.61)-(2.65), one concludes that the solution to Eqs. (2.75)- 
(2.78) must be the image of Eq. (2.39) (i.e., the solution to Eqs. (2.18)-(2.21)) under the mapping Eq. (2.79). 
In other words, the solution to Eqs. (2.75)-(2.78) is 

Uq{j,n) = Q + (v)Uq(j - l,n + 1) + Q-(v)Uq(j + l,n + 1), (j,n) G (2.80) 

By multiplying Eq. (2.80) from left using the matrix U and using Eqs. (2.57) and (2.68), one has Eq. (2.74). 
QED. 
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As a preliminary for the developments in Sec. 3, in the following, important algebraic relations involving 
Q + (^), Q+(^)> an d Q-ty) be extracted from the STI invariance of the a( 4) scheme. 

2.5. Algebraic relations associated with STI invariance 

Let be any given fixed mesh point. Then (j Q + 2 , n G ) G fii and (j Q — 2 , n G ) G fii. Let 

<f(jo + 2,n 0 ), and q(j Q — 2 ,n a ), respectively, be the arbitrary initial data specified at these mesh 
points. Let v 1 ^ 9. Then q(j Q + 1, n 0 + 1) can be uniquely determined in terms of q(j Q , n 0 ) and q(j Q + 2, n 0 ) 
by imposing the basic form Eqs. (2.18)-(2.21) with (j, n) = + 1, n Q + 1), In fact, by using the equivalent 

forward marching form Eq. (2.39), one has 

q(jo + l,n 0 + 1) = Q+(^) q{j 0 + 2 ,n 0 ) + Q_(^) <?(jo, n 0 ) (2.81) 

Similarly, by imposing the basic form Eqs. (2.18)-(2.21) with (j, n) = (j Q — 1, n 0 + 1), one has 

<f(jo - l,n 0 + 1) = Q+(v)q(jo,n 0 ) + Q-(y)q(j 0 - 2,n 0 ) (2.82) 

Among the conditions imposed above, (i) Eqs. (2.18) and (2.19) with (j, n) = (j 0 — l,n G + l) repre- 
sent two conditions linking q(j Q — l,n 0 + l) and q(j 0 -,^o)\ and (ii) Eqs. (2.20) and (2.21) with (j, n) = 
(j 0 + 1, n 0 + 1) represent two conditions linking <f(j 0 + 1, + 1) and g(j 0 , ^o)- Thus the four mesh variables 

in q(j 0 )Ti 0 ) can be determined in terms of q(j Q + 1, n 0 + 1) and <f(j 0 — 1, n Q + 1) by using the four conditions 
specified in the above items (i) and (ii). In fact, by using the equivalent backward marching form Eq. (2.74), 
one has 

q(jo , n Q ) = Q+(v) q(j 0 + 1, n 0 + 1) + Q-W (FO’o - 1, n 0 -h 1) (2.83) 

By substituting Eqs. (2.81) and (2.82) into (2.83), one concludes that 

[Q+WQ-M + Q-{v)Q + {v ) - I]q{j 0 ,n 0 ) 

+ Q + (u)Q + (u)q(j 0 + 2,n 0 ) + Q-(u)Q^(u)q(j 0 -2,n 0 ) =0 

where I is the 4x4 identity matrix and 0 is the 4x1 null column matrix. Because Eq. (2.84) must be valid 
for any choice of q(j 0 + 2, n G ), q(j Q — 2 ,n G ), and <f(jo,n 0 ), the coefficient matrices in front of these column 
matrices must vanish identically. Thus we have 

Q + (v)Q_(v) + Q-(v)Q+(v) = I (2.85) 

Q+ty)Q+(v) = 0 (2.86) 

and 

Q.(v)Q_(v) = 0 (2.87) 

where 0 is the 4x4 null matrix. As an example, one can prove Eq. (2.85) by substituting into Eqs. (2.84) 
each of the following sets of the initial data: (i) q(j Q + 2, n 0 ) = q(j Q — 2, n 0 ) = 0 and q(j 0 -> n 0 ) = (1, 0, 0, 0)*, 
(ii) q{jo + 2 ,n Q ) = q(j 0 - 2 ,n a ) = 0 and q(j 0 ,n 0 ) = (0, 1,0,0)*, (iii) q(j a + 2 ,n 0 ) = cf(j 0 - 2 ,n a ) = 0 and 
q{jo, n a ) = (0, 0, 1, 0)*, and (iv) q(j 0 + 2, n 0 ) = q(j a - 2, n a ) = 0 and q(j 0 , n Q ) = (0, 0, 0, 1)*. Here c* denote 
the transpose of a 1 x 4 matrix c. 

Similarly, by substituting the backward marching relations 

q(j 0 + 1, n a - 1) = Q+( v) q(j 0 + 2, n 0 ) + Q- (v) q(j a , n a ) (2.88) 

and 

q{jo - 1, n a - 1) = Q+(v) q(j 0 , n a ) + Q-{v) q[j 0 ~ 2, n a ) (2.89) 

into the forward marching relation 

q{jo, no) = Q+(i>) q{j 0 + 1, no - 1) + Q_( v) q(j a - 1, n a - 1) (2.90) 
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one has 


(2.91) 


[Q+(u)Q_(v) + Q_(u)Q + (u) - I]q{j 0 ,n 0 ) 

+ Q+WQ+i") Q{jo + 2, n 0 ) + Q-{v)Q _ {u) q(j a - 2, n 0 ) = 0 

Because Eq. (2.91) must be valid for any choice of q(j Q + 2 , n G ), q(j Q — 2 , n G ), and one concludes 

that 

Q+MQ-W + Q- MQ+M = 7 (2-92) 

Q+MO+M = 0 (2.93) 

and 

Q_fy)Q_fy) =0 (2.94) 

By using Eqs. (2.57) and (2.68), it can be shown that: (i) Eq. (2.85) ^ Eq. (2.92) 


Q-{v)UQ-(y) + Q + {v)UQ+{v) = U 


(2.95) 


(ii) Eq. (2.86) ^ Eq. (2.94) 

Q_(v)UQ+{y) = 0 

and (iii) Eq. (2.87) <=> Eq. (2.93) <=$ 

Q + (v)UQ-( v) = 0 


(2.96) 

(2.97) 


2.6. The dual a(4) scheme 

In the above, the a(4) schemes is defined using only the mesh points £ fl\. Independently, it can also 
be defined using only the mesh points £ where 


SI 2 = {(j, n)|j, n = 0, ±1, ±2, ±3, . . . , and (j + n) is an even integer} (2.98) 

For the current ID case where a structured mesh is used, the a( 4) scheme defined over is completely 
decoupled from that defined over Thus, there is no practical reason to carry out computations using the 
two decoupled schemes simultaneously. 

However, to simplify numerical comparisons between the a (4) scheme and the a (3) scheme which is only 
defined over the mesh point set 


n = 9 1 U9 2 = {(j, n) \j, n = 0, ±1, ±2, . . .} (2.99) 

the numerical results to be presented in Sec. 4 are generated using the “dual” a (4) scheme, i.e., the scheme 
formed from the two decoupled a (4) schemes and defined over Q. By the definition of STI invariance, one 
can see that each form of the dual a (4) scheme is also STI invariant. 

3. von Neumann analysis 

Let Gfy, 0) be a 4 x 4 nonsingular complex matrix function of v and the phase angle 0 such that 

q(j , n ) = elj6 [Gfy, 0)] n 6, (j, n) £ Jli; — oo < 9 < +oo; i = \f—\ (3.1) 

is a solution to Eq. (2.39) for all possible complex constant 4x1 column matrices b. (Note: because 
[G(v,Q)] n d = j[Gfy,0)] -1 j for an integer n < 0, [G(v,6)] n is not defined if n < 0 unless [Gfy,0)] _1 exists, 
i.e., Gfy, 0) is nonsingular.) By substituting Eq. (3.1) into Eq. (2.39), one has 

[G{v, 0 ) - e ie Q+(v) - e~ ie Q- (u)] [G( v, 6)] n 6 = 0, n = 0, ±1, ±2, . . . (3.2) 
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Because (i) [G(v, 0)]° = I, and (ii) b can be any complex constant 4x1 column matrix, Eq. (3.2) implies 
that 


G(y,0)=e*Q + {v)+e-*Q-(y) 

By definition, G(v, 0) is the amplification matrix of the forward marching form of the a( 4) scheme. 
By using Eqs. (2.95)-(2.97), one can show easily that 


(3.3) 


U[e i0 Q-(v) + e~ i0 Q+{v)]U [e i0 Q+{v) + e~ i0 Q-(v)] 

= [e i0 Q+(v) +e- i0 Q-(v)\U [e i0 Q-{v) + e~ i0 Q+(u)]U = I 
Thus G(i/, 0) defined in Eq. (3.3) is nonsingular and its inverse is 

[GM )]- 1 = u [e i0 Q-{v) + e~ i0 Q + {v)\ U 


(3.4) 


(3.5) 


Indeed, with the aid of Eq. (2.68), Eq. (3.5) is what one obtains after substituting Eq. (3.1) into the backward 
marching form Eq. (2.74). Moreover, by using Eqs. (2.57), (3.3), and (3.5) along with the fact that Q+(v) 
and Q-(v) are real matrices, one arrives at the important conclusion 

[G( v, 6>)] _1 = t/G(i/,6>)E7 _1 (3.6) 


Hereafter M denotes the complex conjugate of any matrix M. 

For each (i /,0), the four eigenvalues G(v, 0) will be denoted as o^iy, 0), £ = 1, 2,3,4, and referred to as 
the amplification factors of the a( 4) scheme. By using Eq. (3.6), next it will be shown that 


{ 0-1(1/, 6 )' 


a 2 {v,0)' 0-3(1 ',0)’ 0-4 ( 


} = {o"i(M), (?2(v,0), 03 (v,6), o- 4 (i/, 0)j 


(3.7) 


Hereafter z denotes the complex conjugate of any complex number z. 

As a preliminary, first we introduce the following matrix theorems: 

Theorem 1. Let A be a nonsingular N x N matrix with the eigenvalues A^, £ = 1,2 , . . . , TV. Then (i) 
\i ^ 0 , l = 1 , 2 , . . . , IV; and (ii) the eigenvalues of A -1 are 1/A^, l = 1, 2 , . . . , N. 

Theorem 2. Let A be a N x N matrix with the eigenvalues A^, £ = 1, 2, . . . , N. Then the eigenvalues 
of A, the complex conjugate of A, are A^, £ = 1, 2, . . . , N. 

Theorem 3. Let A and B be two similar N x N matrices, i.e., there exists a nonsingular N x N matrix 
S so that B = S~ 1 AS. Then A and B have the same eigenvalues, counting multiplicity. 

The proof of Theorems 1 and 2 is given in Appendix A of [72] while that of Theorem 3 is given on p. 45 of 

[74]. 

To prove Eq. (3.7), consider any (v,0) with v 1 ^ 9. Then because the eigenvalues of the nonsingular 
matrix G(v, 0) are <r^(zv, 0), £ = 1,2, 3, 4, Theorem 1 implies that: (i) a^(v,0) ^ 0, £ = 1,2, 3, 4; and (ii) 
the eigenvalues of [G(v, 0)] _1 are 1/ (Ti(v ,0), £ = 1,2, 3, 4. Next, by using Theorems 2 and 3, and the fact 
that (U^ 1 ) -1 = U, one can see that the eigenvalues of the matrix on the right side of Eq. (3.6) are <r^( v, 0), 
£ = 1, 2, 3, 4. Thus Eq. (3.7) now is an immediate result of Eq. (3.6). QED. 

An immediate result of Eq. (3.7) is 


111 
cri (v,0) <t 2 (M) 0-3(1/, 0) 


1 

0-4(1/, 0) 


= O-!(l/,0) • 0- 2 (l/, 0) ' O- 3 (l/,0 ) ■ 0 - 4 ( 1 /, 0) 


i.e., 


(7l(l/,0)\ ■ \(72(l/,0)\ • \(J3(v,0)\ ■ |ct 4 (M)| = 1 


(3.8) 
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For any given v, stability of the a( 4) scheme requires that 


\oi(v,0) |<1, i= 1, 2, 3, 4 

Thus Eq. (3.8) implies that, for any given v , the a( 4) scheme must be neutrally stable, i.e., 

MM) 1 = 1, ^ = 1,2, 3, 4 


(3.9) 


(3.10) 


if it is stable. As such, Eq. (3.6) does not imply neutral stability of the a (4) scheme. However, it does imply 
that the scheme can only be neutrally stable (i.e., non-dissipative) if it is stable. Here we have reached this 
conclusion without using the explicit form of 0), £ = 1, 2,3, 4. 

4. Numerical results 

To assess the accuracy of the a (4) scheme, consider the model problem with the PDE 


and the exact solution 
We have 


du du _ 
dt + dx 


u = u e (x , t ) = f sin [2tt(x — i)] 


a = X = T= 1 


where A = wavelength and T = period. Let (i) 

/ ,\ def dlle^X^t^ , v def 9 U e {x,f) . * def d 3 u e (x,t) 

r ^xe\P^’)f) — ) 'U'xxe{% ■> 5 ^ n d ^xxxeip^if) — Qx^ 

and (ii) the spatial domain of unit length be divided into K uniform intervals. Thus 

ax = 1 / K, At = vax and t = nAt 


(4.1) 

(4.2) 

(4.3) 

(4.4) 

(4.5) 


where n = number of time steps, and t = total marching time. 

It has been shown numerically that the a( 4) scheme is stable if 

M < 1/3 (4.6) 

On the other hand, (i) the a scheme is stable if 

M < 1 (4.7) 

and (ii) a (3) scheme is stable if 

M < 1/2 (4.8) 

In Tables 1-4, the numerical errors of several computations using the a( 4), a(3), and a schemes are 
presented in terms of the parameters 


E(K, n, v) 


def 


1 


E x (K,n, v) 


def 


i 


1 K ~ 1 
j=0 

(4.9) 

K~ 1 

f ^ [( U b) 7 u xe(%j,t n )] 2 
j=o 

(4.10) 
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(4.11) 


and 


E xx {K,n, v) 


\ 


K - 1 


1 x > 

~J£ ^ ^ [(^xx)j' — u xxe{%ji t n )] 2 


3=0 


E xxx (K,n, v) d = 


N 


K - 1 


1 X X 

^ ^ [(^xxx)j^ ^xxxe(^jj^ n )]^ 


j=0 


(4.12) 


The numerical errors of several simulations with v = 0.1 and t = 9.876 are given in Table 1. For the 
a scheme, as the values of K and n become larger, the values of E and E x are both reduced by a factor 
of about 4 as both K and n double their values, i.e., the scheme is clearly 2nd order in accuracy for both 
u'j and (ux )™ . For the a (3) scheme, the values of I?, E x , and E xx are reduced by factors of about 16, 16, 
and 4, respectively as both K and n double their values. Thus, for this case, the a (3) scheme is 4th order 
in accuracy for both u'j and (u x )™ while only 2nd order in accuracy for (w xx )™ . On the other hand, for the 
a(4) scheme, the values of E , E x , E xx , and E xxx generally are not reduced by constant factors as both K 
and n double their values. However, by comparing their values before and after both K and n increase by a 
factor of 8, one can show that, for the a (4) scheme, the estimated orders of accuracy for u (w x )™, (^xx)j , 
and (w xxx )j are 4.61, 1.83, 1.61, and —0.177, respectively. 

In Table 2, the cases considered have v = 0.1 and t = 10.00 = 10T. For these cases where t is an integer 
multiple of the period T, it is seen that (i) the a scheme is again 2nd order in accuracy for both u'j and 
('U x )™; and (ii) the a (3) scheme is 4th order in accuracy for u'j , (w x )™ , and (w xx )j. On the other hand, for 
the a (4) scheme, the estimated orders of accuracy for u'j , (w x )^ , (w X x)j ? and (w xxx )j are 4.30, 2.17, 2.25 and 
0.164, respectively. 

In Table 3, the cases considered have v = 1/3 and t = 9.6. According to Eq. (4.6), v = 1/3 is right at 
the stability boundary of the a (4) scheme. For these cases, it is seen that, aside from round-off errors, the 
numerical values of u'j and (tt X x)j generated using the a (4) scheme are all identical to their exact solution 
values, respectively. However, one observes that the round-off error for u'j grows linearly with n while that 
for (uxx)™ shows signs of nonlinear growth. 

In Table 4, the cases considered have v = 1/3 and t = 10.00 = 10T, For these cases where (i) the value 
of v is right at the stability boundary of the a (4) scheme and (ii) t is an integer multiple of T, aside from 
round-off errors, the numerical values of u'j , (w x )j , (^xx)j , and (u xxx )^ generated using the a (4) scheme 
are all identical to their exact solution values, respectively. However, because of strong growth of round-off 
errors, these highly accurate results become unsustainable as n increases. 


5. Conclusions and discussions 


A thorough and rigorous discussion of a new high order neutrally stable CESE solver of Eq. (1.1) has 
been presented. Because this two- level explicit scheme is associated with four independent mesh variables 
and four equations per mesh point, it is referred to as the a (4) scheme. As in the case of other similar CESE 
neutrally stable solvers, the a (4) scheme enforces conservation laws locally and globally, and it has the basic, 
forward marching, and backward marching forms. Assuming v 2 ± 9, these forms are equivalent and satisfy 
the STI invariant property defined in Sec. 2. 

Based on the concept of STI invariance, a set of algebraic relations (Eqs. (2.95)-(2.95)) involving the 
coefficient matrices Q + ( v) and Q-{y ) is developed in Sec. 2. As it turns out, these relations can be used to 
construct a simple proof for the fact that the a( 4) scheme is neutrally stable (i.e., non- dissipative) when it 
is stable. Numerically, it has been established that the scheme is stable if the Courant number \v\ < 1/3. 

It is shown in Sec. 4 that the a (4) scheme can be more accurate than the 4th-order non- dissipative a (3) 
scheme, at least for the primary mesh variable u'j. However, the a( 4) scheme has the disadvantage that its 
stability bound is lower than that of the a(3) scheme which is neutrally stable when \v\ < 1/2. 

The CESE development has been driven by a basic idea that each practical scheme be built from a non- 
dissipative core scheme so that the numerical dissipation can be controlled effectively. As such, development 
of the a (4) and a (3) schemes provides a foundation for the development of other more practical high order 
CESE schemes. 
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Figure 1 . — A surface element on the boundary 
S(V) of an arbitrary space-time volume V. 
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Figure 2. — The SEs and CEs. 
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